Classification of mass and nematic orders in bilayer graphene 
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We here classify of all the possible gapped massive and gapless nematic phases, in bilayer graphene. 
Starting from a simple tight binding lattice model, we review the derivation of effective low en- 
ergy Hamiltonian for AB-stacked bilayer graphene. Various discrete and continuous symmetries of 
the non-interacting theory is discussed. Response of the gapped insulating orders for the spinless 
fermions to the quantizing magnetic field is addressed. A precise definition of the quantum anoma- 
lous Hall state is given. We show, even though there are three time-reversal symmetry breaking 
orders for spinless fermions, only one of them corresponds to anomalous Hall state. Upon con- 
structing a particle-hole doubled 16 component Nambu-Dirac spinor, we recognize all the possible 
mass orders, leading to gaps for the quasi-particle excitations near the Dirac points as well as the 
nematic orders, which on the other hand, splits the parabolic dispersion into two anisotropic Dirac 
like conical ones. Transformation of the order parameters under various continuous and/or discrete 
symmetries is reported. We show there are altogether 16 and 12 fermionic bilinears, respectively 
defining all the gapped insulating and superconducting orders. However, physically they correspond 
to only 8 insulating and 4 superconducting phases. Among the gapped superconductors, three are 
spin-singlet, which include uniform s-wave and two spatially inhomogeneous, translational symme- 
try breaking Kekule superconductors. The triplet pairing exhibits an f-wave symmetry around the 
Brillouin zone. Insulating orders in bilayer graphene fall into two broad categories, translational 
symmetry i) preserving layer-polarized, anti-ferromagnet, spin and anomalous Hall phases and ii) 
breaking spin singlet and triplet Kekule-currents. Besides the gapped phases, there are altogether 
56 fermionic bilinears, describing all the nematic orders in bilayer graphene. 32 of them corre- 
sponds to semi-metallic and remaining 24 to superconducting order parameters, which conform to 
16 semimetallic and 8 superconducting states. Some of the nematic superconducting phases break 
the translational symmetry. We name them nematic-Fulde-Farrell-Larkin-Ovchinikov superconduc- 
tors. A brief review over the role of the finite ranged electron-electron interactions, which can 
destabilize the semimetallic bilayer graphene towards various correlated ground states is presented. 

PACS numbers: 71.10.Pm, 71.10.Li, 05.30.Fk, 74.20.Rp, 71.10.-w, 
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I. INTRODUCTION 

Carbon based materials opened a new frontier in con- 
densed matter physics following the successful fabrication 
of single and bilayer graphene.— The low energy excita- 
tions in mono-layer graphene can be described by mass- 
less, chiral Dirac fermions around the six corners of the 
first Brillouin zone. However, only two out of them are 
inequivalenti^ On the other hand, bilayer graphene with 
Bernal stacking, leads to parabolic dispersions around 
those pointsi^ Placed in a weak quantizing magnetic 
field both monolayer and bilayer graphene exhibit quan- 
tized plateaus of Hall conductivity at integer fillings: 
v — ±(4n + 2) in mono-layer and v = ±(4n + 4) in 
bilayer graphene, with n = 0, 1, 2, • • The extra four 
fold degeneracy of all the Landau levels including the ze- 
roth one arises from the valley and the spin degeneracies, 
and the plateaus at v — ±2 (±4) counts the fact that 
only half of the zeroth Landau level is filled in neutral 
mono(bi)layer graphene. The additional two-fold orbital 
degeneracy of the zeroth Landau level in a bilayer system 
is solely due to the quadratic dependence of the energy 
with the momentum^ In the presence of spatially modu- 
lated fields, even though the Landau level structure dis- 
appears, leading to a continuous spectrum, a finite num- 
ber of states always persists at zero energy. Density of 



states at zero energy is twice the flux enclosed by the 
bilayer graphene system, while it is simply proportional 
to the magnetic flux enclosed by monolayer graphenei^ 

In a neutral graphene (either bilayer or monolayer), 
the valence band is completely filled, whereas the con- 
duction band is totally empty, leading to gapless Dirac 
(monolayer) or Dirac-like (bilayer) quasi-particle excita- 
tions around the Dirac points. In pristine monolayer 
graphene, any weak electron-electron interaction is irrel- 
evant due to the vanishing density of states at the charge- 
neutrality point, corresponding to a large domain of at- 
traction for the non-interacting Gaussian fixed point. 9 
On the other hand, in bilayer graphene the density of 
states is finite ( in fact, constant). Consequently, even 
weak electron-electron interactions can be relevant in a 
bilayer graphenei 10 ' 11 Although, in reality such parabolic 
touching of the bands is achieved at the cost of a fine tun- 
ing, where all the hopping amplitudes other than the in- 
plane nearest-neighbor and out of the plane dimer ones, 
are neglected. One therefore still require finite thresh- 
olds of interaction strength to place the system in or- 
dered phases. Nevertheless, the requisite strength of 
the interactions for ordering in bilayer graphene is ex- 
pected to be sufficiently smaller than that in mono-layer 
graphene, as the remote hopping amplitudes are small. 
A similar parabolic band structure can also be realized in 
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checkerboard and Kagome lattices, with particular hop- 
ping structured Otherwise, depending on the relative 
strength of various components of the finite ranged in- 
teraction (repulsive or attractive), fermions in bilayer 
graphene can condense into plethora of insulating, ne- 
matic and superconducting ground states. 

Strength of the Coulomb interaction in a bilayer 
graphene can be tuned quite efficiently by changing the 
gate configuration. For example, interactions are rel- 
atively long-ranged in the presence of a single gate, 
whereas it becomes considerably short-ranged with an 
additional gate (top and bottom gate).— Therefore, 
one can observe different correlated phases in bilayer 
graphene, upon changing number of gates. Some of 
the ordered phases have already been observed experi- 
mentally; with a single gate, the parabolic band struc- 
ture breaks into two anisotropic Dirac like dispersions, 
known as the nematic order}— However with two gates, 
an anti-ferromagnet gap appears to persist even down 
to zero magnetic field. The size of the gap otherwise 
increases monotonically with the strength of the mag- 
netic field. 15 Yet another insulating order can be realized 
by applying a gate voltage between two layers. It pro- 
duces a density imbalance of the charge carriers among 
two layers and that way leads to a gapped spectrum. 
Such a phase is known as layer polarized stated This 
phase has also been predicted to appear from electron- 
electron interactions The above discussion therefore 
raises the following questions: What are the possible 
ordered phases (either fully gapped or nematic) avail- 
able for electrons in bilayer graphene to condense into? 
How do the ordered phases evolve under the influence 
of quantizing magnetic fields? What are the spin triplet 
and/or singlet superconducting state that fermions in bi- 
layer graphene can pair into? Can any nematic (gapless) 
superconducting state be found in bilayer graphene? 

In bilayer graphene the low energy the Dirac-like 
Hamiltonian is eight dimensional, arising from two lay- 
ers, two valleys and two projections of spin (in mono- 
layer graphene the Dirac Hamiltonian is also eight dimen- 
sional, however two layer degrees of freedom is replaced 
by two sublattices). If therefore, the cooper pairs to be 
accommodated, the Nambu particle-hole doubled Hamil- 
tonian is sixteen dimensional. We here show that there 
are altogether 28, 16-dimensional matrices (Ai) that 
anti-commute with the low energy kinetic Hamiltonian, 
yielding non-zero expectation values of the fermionic bi- 
linears ^Ai 1 ^. This number is restricted by the Fermi- 
Dirac statistics, which here translates to a set of alge- 
braic constrains that the matrices (Ai) need to satisfy. 
These orders upon acquiring finite expectation values 
lead to fully gapped phases. Out of 28 bilinears, 16 de- 
fine the insulating orders: layer polarized order (1), anti- 
ferromagnet (3), quantum anomalous Hall (1) and spin 
Hall (3) insulator, odd and even kekule currents in the 
spin singlet (2 x 1) and triplet (2 x 3) channels. The char- 
acterizations 'odd' and 1 even'' of the Kekule orders reflect 
their transformation under the exchange of two Dirac 



points. Numbers in the parentheses indicate the num- 
ber of matrices required to define the corresponding or- 
der parameter. The remaining 12 bi-linears define 3 spin 
singlet and 1 spin triplet superconducting phases. The 
singlet superconducting orders are conventional s-wave 
and two spatially inhomogeneous Kekule superconduc- 
tors. The latter ones break the translational symmetry 
of the lattice into Kekule pattern and offer examples of 
Fulde-Ferrell-Larkin-Ovchinikov type of superconducting 
phases appropriate to the bilayer honeycomb latticej 18 i 19 
The only gapped triplet superconductor exhibits an f- 
wave symmetry, changing sign six times around the Bril- 
louin zone. 

Besides the fully gapped insulating and superconduct- 
ing orders, there are various possible nematic orders 
which upon acquiring finite expectation values split the 
parabolic dispersion into two anisotropic conical ones, lo- 
cated away from the Dirac points. The separation among 
these two points is proportional to the size of the nematic 
orders. The necessary property of any nematic order is 
that it anti- commutes with only one of the matrices ap- 
pearing in the kinetic energy term, while commutes with 
the other. In conjunction with the above definition, a set 
of algebraic constraint restricts the number of nematic 
order parameters in bilayer graphene to 56. 32 of them 
corresponds to nematic semi-metal while remaining 24 to 
the nematic superconductors. However, they corresponds 
to 16 semimetallic and 8 superconducting states. We also 
show that there can be nematic superconductors, which 
break the translational symmetry of the honeycomb bi- 
layer lattice. We name them as nematic- Fulde-Farrel- 
Larkin- Ovchinikov superconductors. 

Organization of rest of the paper is as follows. In the 
next section we discuss the lattice model of hopping and 
the spectrum of free fermions. Sec. Ill is devoted to ar- 
rive at the effective low energy theory in bilayer graphene. 
There we also discous various discrete and/or continuous 
symmetries of the non-interacting system. Behavior of 
various mass orders for spinless fermions in quantizing 
magnetic field is discussed in Sec. IV. A 16-component 
Nambu-Dirac spinor, preserving the spin rotational sym- 
metry is constructed in Sec. V. Various possible insu- 
lating and superconducting masses are shown in Sec. 
VI. Classification of all the nematic orders in bilayer 
graphene is presented in Sec. VII. We discuss the role 
of the electron-electron interactions in Sec. VIII. A brief 
summary of the present discussion is given in Sec. IX. 
In the appendix we present an algebraic derivation of the 
number of the mass orders in bilayer graphene and their 
associated Clifford algebra. 



II. LATTICE MODEL OF FREE FERMIONS 

To set the problem up, let us begin our discussion with 
a simple tight binding description of the free fermions 
(and for simplicity spinless) in a bilayer honeycomb lat- 
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tice. One can decompose it as 



H = + Hi 



(1) 



Here 



i=i.2 

corresponds to intra-layer hopping among the sites of the 
two triangular sub-lattices. The intralayer hopping am- 
plitude is t ^ 2.5 eV— u^{A) is the fermion creation op- 
erator at one of the triangular sublattices generated by 
the linear combination of basis vectors a\ — (y/3,— l)a 
and a 2 = (0, l)a. v(B) is the fermion annihilation oper- 
ator on the other sublattice, then located at B = A + b, 
with the vector b being either b\ = (l/v3, l)a/2, b 2 = 

(l/v3, — l)o/2 or &3 = (— l/\/3,0)o. j = 1, 2 corresponds 
to two layers. Whereas, the inter-layer hopping Hamil- 
tonian takes the following form 



H 



0.0 



H 



0.1 



u 1 - 

■"0,2! 



where 



#0,1 



t±J2 u \(A) u 2 (A) + H.c, 

A 

t^sJ^Ai^MA-b^ + H., 



A A 



and 



K* = t 1 B 2 B J2 v l(B)v2(B-b l )+H.c. 



(3) 

(4) 
(5) 

(6) 



energy spectrum is composed of four band with disper- 
sions 



Ei(k) = ±- 



1 



t 2 \mi 



(10) 



and 



E 2 (k) = ±-^4t*|/(/0| a + *i 



o(|JWI 4 )- 



(11) 

Note that /(fc) is zero at K = (l,l/y/3)2ir/ay/3 and -K, 
located at the two inequivalent corners of the Brillouin 
zone. Near the Dirac points Ei(k) vanishes and the dis- 
persion is comprised of two parabolic bands, touching 
each other. On the other hand, spectrum of E 2 (k) is 
gapped near the Dirac points, with band gap ~ 2 i^. 4 
Such parabolic degeneracy can only be achieved after set- 
ting the remote hopping terms to zero. Inclusion of re- 
mote hoppings t B 2 B , splits the parabolic bands into four 
Dirac cones, known as trigonal wraping. In general the 
eigenfunctions of E 2 have finite overlap on both the sub- 
lattices. However, in the vicinity of the Dirac points they 
can be considered to be localized on the sites Aj, with 
i = 1,2. Expanding f(k) near ±K yields 



f(±K 



■<f) = 2 ( ±Qx + 1 q y> ■ 



Therefore 



E x (±K + q) 



1 

2ro 



ill 



(12) 



(13) 



with m — 4t±/3t 2 a 2 being the mass of the parabolic 
dispersion, a is the lattice spacing. In bilayer graphene 
m f=a 0.03m P , where m e is the mass of free electrons. 



Various components of the inter-layer hoppings have been 
measured experimentally^ Their estimated strengths 
are t± ~ 0.5 eV, t\ 2 B ~ 0.1 eV, t 1 ^ ~ 0.3 eV. To di- 
agonalize Hq, let us define a 4-component spinor 



(7) 



*fe = \ui(k),vi(k),U2(k),V2(k) 
The tight binding Hamiltonian in this basis is 



H 



where 



tf(k) t± t 12 B f(k)\ 

tf*{k) t l2 B f(k) 

t±. tf*(k) 

t A 2 B f*{k) t l2 B f*{k) tf(k) / 



f( k ) = exp (i k ■ b^j , 

z=l,2,3 



(8) 



(9) 



and f*{k) is the complex conjugate. 

Next we drop the remote hopping amplitudes, t\ 2 B and 
t B 2 B , as they are much smaller than t and t±. Then the 



III. EFFECTIVE LOW ENERGY THEORY 

A. Lagrangian 

In the previous section we show that near the Dirac 
points (as q x , q y — > 0), states in the gapped band (E 2 (k)) 
are localized on sites of A-sub- lattices, A\ and A 2 , named 
as 'dimer sites'. Hence, in the low energy limit these sites 
are of no dynamical importance. In this section we will 
derive the form of the effective low energy Hamiltonian 
after integrating out the high energy band, E 2 (k). Our 
derivation mimics the one performed by O. Vafek in Ref. 
(22). 22 Nevertheless, it is worth reviewing that deriva- 
tion briefly to facilitate further discussion. The parti- 
tion function describing the free motion of the spinless 
fcrmions is 



"J dr (®t+Haa) 1pA+il A HABlpB+1pB H BA*pA 



(14) 
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where ipx = (ipi.x, 1 p2,x) with X = A, B are two com- 
ponent spinors and i — 1,2 correspond to the layer in- 
dices. For convenience we neglect the remote hopping 
amplitudes and keep the one among the 'dimer sites' (t±) 
and the intra-layer hopping (t). In this representation 
Hab, Hbb and Haa are two dimensional matrices and 
read as 



A(~ t\j\t ~ 200 meV) is the high energy cut-off repre- 
senting the range of energies over which the dispersion is 
approximately parabolic. The mutually anti-commuting 
four component Hermitian gamma matrices belong to the 
representation^ 3 - 



H BB 








H AA 



t ± 
t ± 



and 



tf(k) \ _ t 

tf*(k) I ' AB' 



(15) 



(16) 



Even though, Hbb is trivial at the bare level, once we 
integrate out the high energy modes, it gets renormalized. 
The single praticle Green's function for ipA held is 

lUJnh + t±(T x 



G AA {iUn) = 2 2 

n T 1 l 



(17) 



where uj n = (2n + l)irT are the fermionic Matsubara 
frequencies, and T is the temperature. Expanding the 
action to the quadratic order and integrating out the ipA 
field, yields a renormalized partition function for the ?pB 
field 

— $ y^ t[ . i>B( iuJ n)(- iuJ n+ H BB-HABG A A(iu n )H BA )lp B (iu n ) 



(18) 



Since we are interested in the modes near the charge- 
neutrality points, one can simply set uj n = in the last 
term of the effective action. 



B. Hamiltonian 

The imaginary-time non-interacting Lagrangian (Lq) 
is related to the single-particle Hamiltonian (Ho) via 
Lq = —iu n + Hq. Therefore the low energy Hamiltonian, 
quadratic in momentum, is 



Hn = 







T-(q x + iq v ) 2 



T-{<lx-iqy) 2 







(19) 



near one of the Dirac point at K. Taking into account the 
Fourier modes near the other Dirac point at —K, the four 
dimensional Hamiltonian describing the low excitations 
reads as 



H = 72 



2 m 



-2 q x q y 
2 m 



(20) 



in the basis of four-component spinor ^(x), defined as 



d q 



v{(K + q), 



(27TO) 2 

v\{K + q),v\ (-K + q),vl (-K + q) 



(21) 



7o 



<J Z 

G z 



7i 



(Jy 

-a,. 



72 



a x 
a x 



The two remaining matrices can then be chosen as 
ov 



73 







75 = 





i<r„ 







(22) 



(23) 



These matrices satisfy the anti-commuting Clifford alge- 
bra {7^,71/} = 25 Mt/ , for /i, v = 0,1,2,3,5. The form of 
the free Hamiltonian is identical to the one for monolayer 
graphene in quadratic orderi 24 ' 25 The underlying reason 
is as follows. After integrating out the high energy bands, 
the remaining lattice points on B-sublattices also consti- 
tute a honeycomb lattice, preserving the C% v symmetry 
around each site. Therefore, the free Hamiltonian needs 
to be invariant under a rotation by 27r/3 around the Dirac 
points, which restricts the kinetic energy Hamiltonian to 
the announced form. 



C. Symmetries 

The free Hamiltonian respects an emergent chi- 
ral SU C (2) symmetry, generated by {«7o73, «7o75, 73s}, 
where 735 = 17375. The third entity of the group, is 
the generator of translation ! 22 ! 23 A similar chiral sym- 
metry is also present in the emergent low energy the- 
ory of the massless Dirac fermions in graphene^ 3 - and 
d-wave superconductors^ In addition to the chiral sym- 
metry, Hq is also invariant under the exchange of the 
layer indices, as well as the Dirac points. 22 These two 
reflection symmetries are generated by I12 = 72 and 
Ik = 271751 accompanied by the inversions of the mo- 
mentum axis q y — > — q y and q x — > — q x , respectively. 
Besides these, Hq is also invariant under the time re- 
versal symmetry, since it describes motion of fermions 
on a lattice. Time reversal symmetry is represented by 
an anti- unitary operator I t = UK, where U is unitary 
and K is the complex conjugation. In our representa- 



tion U = «7i75 



Ik- Therefore if 



-1, 



as it should be, since It is the time reversal operator 
for the spinless fermions^ Oneself can arrive at the 
same effective low energy Hamiltonian by using the K • p 
approach. 22 ' 27 



IV. SPINLESS FERMIONS AND ANOMALOUS 
HALL STATE 

Before considering spinfull, it is worth understand- 
ing the possible gapped state of spinless fermions in bi- 
layer graphene and their behavior in quantizing mag- 
netic fields. In monolayer graphene, the linear dispersion 
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makes all the short-ranged electron-electron interactions 
irrelevant near the noninteracting Gaussian fixed point. 
The long range Coulomb interaction (~ 1/r) is also ir- 
relevant, but only marginally 23 ' 29 ' 30 On the other hand, 
due to the quadratic band structure, all the short-ranged 
interactions are marginal in bilayer graphene. Hence, 
the semimetal-insulator transitions can take place even 
for weak interactions ! 11 ' 22 Any order parameter, leading 
to a gap in the spectrum, must anticommute with the 
entire free Hamiltonian, so that all the terms enter as 
sum of the squares in the expression for energy. A fi- 
nite gap then exists everywhere in the Brillouin zone. 
For spinless fermions, there are four such candidates: 
(70 , 73 j 75 j £7175), in a bilayer graphene. The first 
three members break the chiral SU C (2) symmetry of the 
free theory down to U c (l), whereas the last one lacks 
the time reversal symmetry, but preserves chiral sym- 
metry. (5 , ^7o^ / ) is the order parameter associated with 
the layer polarized state, leading to an imbalance of the 
electronic density among two layers. i^i^2^} corre- 
sponds to Haldanes circulating current among the sites 
on same sub-lattice. 3 - Otherwise, it propagates in the 
same directions in two layers. The remaining two enti- 
ties, (73,75) vp) break the chiral as well as the time 
reversal symmetry^ One can however, define an anti- 
unitary operator as in Ref. (32), I t — 17375 K, under 
which three chiral symmetry breaking masses are even, 
whereas the original time reversal odd mass, still remains 
odd. Otherwise, if = — 1, therefore does not corresponds 
to true time reversal operator. (73, 75) \&) correspond 
to Kekule currents, which additionally break the transla- 
tional symmetry of the lattice into Kekule pattern. Oth- 
erwise the first (second) member is even (odd) under the 
exchange of two Dirac points. In monolayer graphene 
these two matrices are replaced by 17073 and £7075, giv- 
ing rise to different realizations of the spatially modu- 
lated Kekule bond density waves £2- However, they are 
time reversal symmetric. Next we argue, even though the 
kekule currents break the time reversal symmetry, only 
the Haldane order (47172) corresponds to the anomalous 
Hall insulators. 

Placed in a weak quantizing magnetic field, bilayer 
graphene exhibits plateaus in Hall conductivity at fill- 
ings v = ±4(77 + 1). The orbital effect of the mag- 
netic field can be captured via a minimal substitution 
Qi Qi — M in H in Eq. (|2"D|) . giving i?o[-<4]- Magnetic 
field then reads as B = CijdiAj. The spectrum of l?o[-<4] 
is comprised of a set of macroscopically degenerate Lan- 
dau levels at well separated energies E n = \Jn(n — \)B 2 
with n = 0, 1, 2, • • •. The zeroth Landau level in bilayer 
graphene (with n — 0, 1) carries additional two fold or- 
bital degeneracy due to the parabolic dispersion in the 
vicinity of the Dirac points. States in the zeroth Landau 
level near two valleys (K and —K) reside on the comple- 
mentary layers, 1 and 2 respectively, for example. These 
two sets of zero energy state constitute a two dimensional 
basis "Ho- Any matrix that commutes or anticommutes 
with the Hamiltonian Hq[A], leaves the space invariant. 



There are four matrices, falling into the second category, 
{70, 73, 75, *7i72}- Together they also close a Cl(3) X U(l) 
algebra of the order parameters, where the U(l) part is 
constituted by the last entry^ To understand the be- 
havior of these orders in the presence of magnetic fields, 
let us consider an auxiliary Hamiltonian 

H [mi, 7712, 7773, 777,4] = 

H [A] + TO170 + TO273 + TO375 + m 4 i7i7 2 . (24) 

The eigenvalues of H[mi, 777,2, 1713, 0] are at 
±\/n(n — 1)£> 2 + to 2 + m| + m§. Any linear com- 
bination of mi, 777,2, TO3 reduces the chiral SU C (2) 
symmetry to a U c (l). Therefore in the presence of the 
chiral symmetry breaking orders the zeroth Landau level 
splits to Eq = ± \/m\ + m 2 , + m 2 , whereas Landau levels 
at finite energies are only shifted. Hence it is always 
energetically advantageous for the system to develop 
such mass orders in the presence of the magnetic fields, 
to maximally lower the energy. The mechanism of devel- 
oping a chiral symmetry breaking mass order is known as 
magnetic catalysis, discussed previously in the context of 
Dirac fermions subject to magnetic fieldsj 35 ' 37 It can also 
be the underlying mechanism behind the formation of 
Hall states at fillings v — 0, ±1 in monolayer graphene. 36 
However the energy spectrum of H[0, 0,0, 7774] within 
the zeroth Landau level is sign(m4)m4. Rest of the 
Landau levels are shifted to ±y/n(n — 1)B 2 + m\. A 
finite 7774 therefore shifts the zeroth Landau level, can 
only be realized by changing the filling factor from 

the neutrality (v — 0). Concomitantly, it leads to 

2 

quantized Hall conductivity a xy — ±2^-. On the other 
hand, 777,2 and 7773 even though break the time-reversal 
symmetry, system discerns zero Hall conductivity when 
7772 7^ and/or 777,3 7^ 0. Therefore, only the Haldane 
order (77172) corresponds to quantum anomalous Hall 
insulator. A similar conclusion can be arrived at upon 
computing the expectation values of aforementioned 
orders, with the non-interacting wavefunctions of the 
zeroth Landau leveL 38 It can also be justified from 
that fact that M c y = (70 , 73 , 75) transform as a 
vector under the chiral SU C (2) rotation generated by 
{«7o73, £7o75j*7375}, whereas M c . s = 77172 is a scalar 
under under the chiral transformation. Hence all the 
three chiral symmetry breaking orders must lead to 
identical spectrum of Hall conductivity in quantizing 
magnetic fields, v — Hall state. 

Recently, there has been various proposal for the elec- 
tronic ground state in bilayer graphene subject to mag- 
netic fields In particular, a recent experiment shows 
that bilayer graphene exists in a gapped phase in the ab- 
sence of magnetic field when the system is dual gatedJ^ 
Otherwise the gap grows monotonically the magnetic 
field, whereas it gradually disappears with a weak per- 
pendicular electric field. These observation predicts that 
the zero field order is likely to be an anti-ferromagnet 
and rules out the possibility of a layer polarized state. 
However, for week fields the most promising candidate 
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of the correlated state in that system is a partially 
spin polarized state, with coexisting easy plane anti- 
ferromagnet and and easy axis (in the direction of the 
field) magnetization^ 



V. NAMBU-DIRAC FERMION 

In the previous section, we considered the insulating 
orders for spinless fermions in a bilayer graphene sys- 
tem, leading to gaps in the spectrum. We name them 
as mass orders. We study their response to quantiz- 
ing magnetic fields and modification of the Landau level 
spectrum when those orders acquire an non- vanishing ex- 
pectation values. Next we wish to find all the mass or- 
ders, including the insulators as well as the superconduc- 
tors and various nematic orders. To accommodate all 
the order parameters, we need to introduce a particle- 
hole doubled 16-componcnt Nambu-Dirac fermion, de- 
fined as * = (* p ,*, l ) T > with * P = (* P ,t,*p,i) T and 
* P = (*h : ;,-*/ l , t ) T with 



v 1<(T (K + q),v 2 ,a(K + q),Vi t(T (-K + q),v 2 , a {-K + q) 



generated by h 2 = r <g) er ® 72 and I K = t <g) er ® «7i75 ■ 
The generator of translation is P = r 3 (g> cr <g 27375 • Ho 
also commutes with the number operator N = T3<g<7o<gi4. 
One advantage of this representation of the Nambu-Dirac 
spinor is that three generators of the rotation of electrons 
spin assume a simple form S = tq <g a <g 1^. The free 
Hamiltonian {Ho) commutes with S. 

Besides the above 16-component Nambu-Dirac spinor, 
we also define an unitarily equivalent one $ = ($ p , $^) T , 
where $ p = ($ P ,t, *W|,) T and $ h = ($ h , t , $/ l , i ) T , with 



HM = 



(K - q),v{jK - q),vU-K - Q)Aj~K - 



similar to a Nambu-Dirac spinor recently considered in 
the context of mono-layer graphene4i In this basis, the 
tight binding Hamiltonian in the low energy approxima- 
tion then takes the form 



(27) 



where Hq in the first quantization reads as 



Ho 



T 3 ® <J 



72 



2 m 



7i 



-2 q x q y 
2 m 



(28) 

where the four dimensional gamma matrices belong to 
the aforementioned representation. The effective mass 
of the quasiparticle excitations is m = t±/v F , where 
vp = \/3 ta/2 is the Fermi velocity in single layer 
graphene. The two component Pauli matrices (to,t) op- 
erates on Nambu's, whereas (o~o,&) on the spin indices. 
Before we proceed to identify the massive and nematic 
order parameters, it is worth pausing to register the sym- 
metries of the quadratic Hamiltonian H Q . The reflection 
symmetries of Hq under the exchange of the layers and 
the Dirac points, mentioned in Sec. II, are respectively 



(29) 



and 



v\JK - q),vl a {K - q), vlj-K - q), v\ M {-K - q) 

(30) 



This representation ($) is unitarily related to the previ- 
ous one ( 1 J r ) 



Ho) 



(25) 

Q) 
(26) 



h © (io- 2 <g 72) 



(31) 



where ig is a 8-dimensional unity matrix. Soon we will 
appreciate the usefulness of such unitary equivalence. A 
similar definition of spinor has recent been used to de- 
scribe all the possible masses in monolayer graphene. 42 
We can immediately notice that there are altogether 64 
matrices (Ai), where 



M = (t ,t 3 ) (gi (<j , a) ® (7o,73,75,i7i72) , (32) 



and 



M = (ti,t 2 ) <g ((To, a) <g (^4,«7o73,«7o75,«7375) , (33) 

which anti-commute with the kinetic energy Hamiltonian 
Hq in Eq. (f28|) . Naively, one may therefore assume that 
all the 64 fermionic bilinears of the form ^AA^ will lead 
to gapped quasi-particle excitation. However, such abun- 
dance is clearly an artifact of Nambu doubling of original 
degrees of freedoms. Next we show that this number is 
drastically reduced by some algebraic constrains that the 
mass order parameters need to satisfy. Upon imposing 
the constraints, we show that only 28 bilinears lead to 
gap in the excitation spectrum. 



VI. MASS ORDERS 

Next we wish to derive the algebraic constraints that 
all the mass order parameters need to satisfy. We derive 
them separately for insulating and the superconducting 
orders. 
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A. Gapped insulators 



7. (T3 ® ® 17172) 1 S r ) : spin Hall insulator 



All the insulating order parameters commute with the 
number operator (N). Therefore, one can have either 
tq or T3 in the Nambu space. Otherwise, in general any 
insulating order (INS) is restricted to the following form 



INS 



Mi 


' 





M 2 _ 



(34) 



where Mi and M2 are 8-dimensional Hermitian matri- 
ces. Let us first consider the insulating orders with To in 



INS 




M 2 


and 


Mi 








Mi 



h 








-iR 



i? x Mi x i? $ ft 



& (Mi - (i? x Mi x 



where i? = 02 ® 72- In our representation, all the ma- 
trices are either purely real or purely imaginary. Hence 
either M, 7 = +Mi or Mj = -Mi, respectively. To- 
gether with this condition, (02 <8> 72 ) T = — 02 <8> 72 gives 
us the requisite condition for non-zero expectation value 
of the gapped insulating orders ((INS) ^ 0), 



(02 53 72) x Mj x (cr 2 (Xi 72) = ~Mi. 



(35) 



Only the following bilinears (and hence mass order 
parameters) meet the above criteria, 

1. (t <x> cr ® 7o) V I / ) : layer polarized state, 

2. ('ft (tq <E) (To ® »7i7 2 ) 1 &) : anomalous Hall insula- 
tor, 

3. (to ® cto <8) 75) V I / ) : odd spin singlet Kekule cur- 
rent , 

4. (\E ,t (to <X> a ® 73) ^) : even spin-triplet Kekule 
current. 

The classification even, and odd reflects the transforma- 
tion of the order parameter under the exchange of Dirac 
points. On the other hand, with T3 in the Nambu space, 
i.e. Mi = — M2 the above condition in Eq. (|35|) reverts 
to 

(cr 2 <8> 72) x M7 x (ct 2 ® 72) = +Mi. (36) 
It allows us to capture the remaining insulating masses: 



5. (T3 (8 (To ® 73) ^) : even spin-singlet Kekule cur- 
rent, 

6. (T3 (g> <r €5 70) ^) : antiferromagnet, 



8. (^ (T3 ® <r ® 75) ^) : odd spin-triplet Kekule cur- 
rent. 

Hence, there are 8-insulating orders, towards which the 
semimetallic bilayer graphene may suffer instabilities and 
open gap at the Dirac points. As one can see that it re- 
quires 16 linearly independent matrices, anti-commuting 
with the kinetic energy Hamiltonian, to define all the 
gapped insulating orders. However, some of them are 
connected by the ordinary rotations of electrons spin. 
All the spin-singlet orders, e.g. layer polarized, quan- 
tum anomalous, singlet Kekule currents break the dis- 
crete layer- inversion symmetry, generated by /12. On 
the other hand, the spin triplet orders, anti-ferromagnet, 
spin-Hall, spin Kekule orders additionally break the spin 
rotational symmetry. Thus, such ordered phases are al- 
ways accompanied by 2 massless Goldstone modes. In 
the low energy limit, there is an internal U(l) symme- 
try among various realizations of Kekule orders (for both 
spin singlet and triplet), (TO173, 771275). However, such 
an emergent internal symmetry is clearly an artifact of 
the parabolic band approximation of the dispersion near 
the Dirac points. The actual lattice always reduces such 
an emergent symmetry. Hence in an ordered phase with 
finite Kekule current there is no massless mode arising 
from the breaking of the internal U(l) symmetry. 



B. Gapped superconductors 

Next we find the gapped superconducting states avail- 
able for fermions in bilayer graphene to pair into. Super- 
conducting or the off-diagonal order parameters, on the 
other hand anticommute with the number operators (JV) 
and takes the following form 



SC 






M 


Mt 






(37) 



in general. M is an 8-dimensional matrix. SC can also 
be casted in the following form 



SC 







M 


Mt 






h 








-iR 






-%M R 


iRM^ 






(38) 



Therefore, any off-diagonal order with nontrivial expec- 
tation value must satisfy the following condition 



(cr 2 ® 72) M T = M (cr 2 ® 72) ■ 



(39) 



The superconducting orders concur with this constrain 
are 



1. (ty [(ti cos (f> + t 2 sin ( 
singlet s-wave, 



c ® 17073] *> : 



spin- 
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2. (\& [(n cos + r 2 sin 0) ® a (8 17075] : spin- 
triplet f-wave, 

3. (\E ( [(n cos + r 2 sin 0) ® oo ® ^4] : spin-singlet 
s-Kckule, 

4. (\fr [(n cos + r 2 sin 0) ® cr ® *737s] : spin- 
singlet p-Kekule, 

superconductors. Here <p is the superconducting phase. 
The s-wave order is even under the exchange of layers 
and Dirac points and translationally invariant. Similar 
s-wave phase can also be realized in neutral mono-layer 
graphene if the onsite attractive interaction is sufficiently 
strong* 4 ^ The translationally invariant f-wave SC order, 
is odd (even) under Dirac point (layer) exchange and 
changes its sign six times around the Brillouin zone, sim- 
ilar to the one in mono-layer graphene. 44 The Kekule su- 
perconductors, on the other hand, break the translational 
symmetry of the lattice into Kekule pattern, odd under 
the exchange of layers and spin singlet. Otherwise, s(p)- 
Kekule is even (odd) under the Dirac point exchange. In 
a monolayer graphene, the Kekule superconductors are 
spin-triplet4i 

Therefore, fermions in bilayer graphene can pair into 4 
gapped superconducting phases. However, one requires 
12 linearly independent matrices to define all of them. 
Hence, altogether there are 

16 (insulators) + 12 (superconductors) = 28 

fermionic bi-linears which can give rise to gap in the 
quasiparticle dispersion, near the parabolic degeneracy 
points. 

VII. NEMATIC ORDERS 

Besides the fully gapped gapped orders, fermions in bi- 
layer graphene can also develop interaction induced ne- 
matic orders. Before we proceed to recognize all possible 
nematic orders, it is worth pausing to provide a general- 
ized definition of the nematic orders. 



A. Definition 

For simplicity let us consider spinless fermions once 
again. Recall that in the vicinity of the parabolic touch- 
ing points the kinetic energy Hamiltonian takes the form 

tfofe, ?1 ,]=7 2 (^) +71 (40) 
\ 2m I \ 2m J 

A nematic order is associated with a matrix (M. ) which 
anti-commutes with only one of the matrices appearing 
in the free Hamiltonian (namely, 71 and 72 ), while it 
commutes with the other one. If such an order develops 
a finite expectation value then the parabolic bands split 



into two Dirac cones separated by an amount propor- 
tional to the size of the nematic order. Otherwise, the 
nematic order is invariant under rotation by tt. Note, the 
parabolic bands are associated with a Berry phase 27T. 
Therefore, an order parameter can in principle split the 
parabolic band into two Dirac cones, each of them carry- 
ing a Berry phase n. To gain further intuition about the 
nematic orders let us consider two specific examples: (1) 
M = A171 and (2) M = A 2 7 2 . Here Ai i2 correspond 
to sizes of two nematic orders. In the former situation, 
the Dirac cones appear at q x = q y = ±v / mA7/2m, 
whereas a finite A 2 , gives birth to two Dirac cones at 
q x = 0,q y = ±V2A 2 m/2m. 

B. Nematic semimetals 

Next we consider the 16-component Nambu-Dirac 
spinor (vp) and wish to write down all the possible ne- 
matic semi-metals (commutes with the number oper- 
ator) as well as superconductors (anti-commutes with 
the number operator). Natively, one can assume there 
are all together 64 semimetallic nematic order parame- 
ters, of the form (t ,t 3 ) ® (cro,ff) ® Nm, where Nm € 
(7i>72, i7o7i)*7o72, «7i73, «7273, «7i75, 27275)- However, 
to acquire finite expectation values, the semi-metallic 
nematic orders need to satisfy one of the constrains in 
Eq. ([35]) (if spin-singlet) or Eq. {36) (if spin-triplet). 
In conjunction with the above definition, those con- 
strains restrict the possible spin-singlet nematic insula- 
tors ((^M*) ^ 0) to the following bilinears: 

1. (* t {T ®a (8)71, t ®ct ® 72)*), 

2. (t 1 * (r 3 ® <r ® 17072, f3 ® o- ® «7o7i) *)i 

3. (* t (t <E> (To ® »7273, t ® CT ® «7l73) 

4. (*t (r 3 cr ® 17375, T3 ® (Tq ® 17175) 

Elements within the same group are connected by ro- 
tation of 7r/2 around the Dirac points, generated by 
Irk = To ® o~o ® *7i72- The first member of each group 
is odd under exchanges of two layers, whereas the second 
entry is even under the same operation. This symmetry 
is generated by I12 = t ® 00 <£> 7 2 . The nematic orders in 
last two group (3 and 4) break the translational symme- 
try generated by P — r 3 ® oo ® 17375, where those in the 
first two groups (1 and 2) are invariant under translation. 
Under the exchange of two Dirac points, the first member 
of the group 1 , and the second members of the remaining 
groups, are even. Whereas the other members of each 
group are odd under that interchange. One can write 
down all the bilinears (order parameters) corresponds to 
the spin-nematic orders, 

1. (* f (r 3 ® a ® 71, r 3 ® a ® 72)*), 

2. (*t (tq®B® ijol2, to® a® 17071) 

3. (* t (r 3 ® a ® «7 2 73,T3 ® <?® 17173) 
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4. (\E rt (t q a (g 17275, To <g <? Hi 75) 

Transformations of the spin-nematic orders under the ex- 
change of two layers, Dirac points, translational are iden- 
tical to the ones for spin-singlet nematic orders, since all 
the symmetry operator have cto (two dimensional iden- 
tity matrix) in the spin block. Apart from the massless 
Dirac fermionic excitations, spin-nematic phases are also 
accompanied by two massless Goldstone modes in the 
ordered states. 



C. Nematic superconductors 

Other than the nematic insulators fermions can also 
pair into various nematic superconductors. Though, one 
can find 64, Nambu-Dirac bilinears, which anti-commutes 
with the number operator and only one of the matrices 
in the free Hamiltonian, only the following pairing order 
parameters satisfy Eq. (|39l) 

1. (*t [(n cos + t 2 sin (ft) <g er <g 71] *> 

2. (<Jrt [(ti cos 4> + t 2 sin (ft) <g a <g 72] *) 

3. ("ft [(ti cos <j) + t~2 sin (ft) <g ao (g 17071] 

4. [(ti cos + r 2 sin (ft) (g (To 17072] *) 

5. (^t [(ti cos + r 2 sin (ft) <g <? <g 17173] 

6. [(ti cos (j) + t 2 sin 0) <g <? (g ^7273] 

7. (\E rt [(n cos + r 2 sin </>) !g ao ig 17175] 

8. [(ti cos0 + r 2 sin0) !g cr <g> «727s] 

where ^ is the superconducting phase. Last 4 supercon- 
ducting orders preserve the translational symmetry and 
the cooper pairs are formed by pairing fermions with 
momentum K + q and — K — q or q — > —q. The re- 
maining 4 parings are spatially inhomogeneous with pe- 
riodicity 2K and breaks the translational symmetry of 
bilayer honeycomb lattice into Kekule pattern. Cooper 
pairs in those channels are formed by gluing fermions 
with momenta K + q and K — q and K — > —K. We 
name them as nematic-Fulde-Ferrell-Larkin-Ovchinikov 
superconductors. The 1st, 4th, 6th, 8th nematic super- 
conducting pairings are even, while the remaining four 
pairing are odd, under the exchange of two layers. On the 
other hand, the 1st, 3rd, 5th and 8th superconducting or- 
ders changes sign under the exchange of two Dirac points, 
while the rest four pairings remains invariant. Elements 
of the pairs (l st , 2 nd ) , (3 rd , 4 th ) , (5 th , 6 th ) , (7 th , 8 th ) of 
the nematic superconducting orders transforms into each 
other under the tt/2 rotation around the Dirac points, 
generated by Irk- 

Above we present all the possible nematic, including 
semi-matallic as well as superconducting orders. All to- 
gether there are 

32 (semimetals) + 24 (superconductors) =56 



fermionic bi-linears that define all the nematic orders in 
honeycomb bilayer graphene. Note, one can write a set 
of matrices as Nem=iHlM.k, where Hq is one of the 
two matrices appearing in the kinetic energy Hamilto- 
nian, whereas Mk is one of the 28 matrices, defining 
the massive order parameters. By construction, Nem is 
Hermitian, anti-commutes with one of the matrices in 
free Hamiltonian, while commutes with the other one, 
hence meets the definition of the nematic order parame- 
ters. Therefore, the total number of nematic semi-metal 
is 2 x 16 (number of gapped insulators) = 32, and that of 
superconductors is 2 x 12 (number of gapped supercon- 
ductors) = 24, in accordance with our explicit computa- 
tion, yielding total 56 nematic order parameters. How- 
ever, physically they corresponds to 16 semimetals and 8 
superconductors . 



VIII. INTERACTIONS 

In this section, we present a qualitative discussion 
on the role of electron-electron interactions in bilayer 
graphene. As mentioned previously, all the four fermion 
interactions are marginal in the bare level. That al- 
lows one to perform an weak coupling expansion about 
the symmetric semi-metallic ground state to study its 
instabilities towards the formation of various correlated 
phases. The interacting theory in bilayer graphene can 
be expressed in terms of 18 quartic interactions, 22 How- 
ever, not all the 18 coupling constants are linearly in- 
dependent. There exist a set of linear constrains, so 
called Fierz identity, which allows one to write each of the 
quartic terms as a linear combination of the others . 23 > 45 
Such linear constrains restrict the number of indepen- 
dent quartic terms to nine (9). For example, one can 
write all the interactions in the spin-triplet channel as 
linear combinations of the singlet ones. 

With repulsive Hubbard interaction (U), a one loop 
renormalization group calculation shows that system 
finds itself in a state with a staggered pattern of spin 
among the two layers. This prediction can also be jus- 
tified from the strong coupling physics. 22 It is expected 
that each of the layers is anti-ferromagnetic ordered at 
least when t± = and U/t 3> l. 46 However, the rela- 
tive orientation of the anti-ferromagnet order in two lay- 
ers, is arbitrary, when t±_ = 0. Upon turning on tj_, 
the sub-lattice magnetization on two layers assumes a 
staggered pattern. Using a similar argument, one can 
also predict the possible ground state if the repulsion 
(V2) among the fermions living on the next-neighbor sites 
on same layer is the strongest component of the finite- 
ranged Coulomb interaction. If the layers are completely 
decoupled (t± = 0), each layer is expected to find it- 
self in the quantum spin Hall insulator phase, at least 
when Vijt » l4£ This phase supports circulating cur- 
rents among the sites of the same sub-lattice. Otherwise 
it orients in opposite direction on two sublattices^ Its 
orientation is opposite for two spin projections. The spin 
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Hall insulator preserves the total time reversal symme- 
try as well as the inversion symmetry. The spin Hall 
insulator additionally breaks the spin rotational symme- 
try. Therefore the ordered phase is accompanied by one 
massive and two massless modes ; 47 ' 48 When t± = 0, the 
orientation of the Haldane's circulating current in two 
layers are completely independent. A small t±, however 
locks the circulation in two layers in the same direction. 

When the interaction is relatively long ranged possi- 
bly an unconventional phase, nematic order arises. Un- 
like the fully gapped phases, a nematic order splits the 
parabolic band into two Dirac cones. However they ap- 
pear at different location in the Brillouin zone than the 
Dirac points, which on the other hand remain gapped. 
The separation among these two cones are proportional 
to the magnitude of the nematic order. The experimen- 
tally observed nematic phase appears to be translational 
symmetric^ Note that any lattice model with density- 
density interaction contains both intra-valley (forward) 
as well as inter-valley (backward) scatterings. Their rel- 
ative strength however depends on the range of the in- 
teraction. Therefore one can find a rich phase digram of 
various correlated phases simply by tuning the relative 
strength of these two scatterings 

If on the other hand, the net interaction acquires an 
attractive component, fermions in bilayer graphene may 
condense into variety of superconducting states. An at- 
tractive interaction can arise, for example from electron- 
phonon interactions or novel proximity effect. An on- 
site attraction can favor a spin singlet s-wave supercon- 
ducting order, as in graphenei 43 ' 49 Attractive interac- 
tion among the fermions living on the same layer, but 
at the next neighbor sites can support a spin triplet f- 
wave superconducting state. The superconducting order 
parameter changes its sign six times around the Brillouin 
zone, similar to the one appropriate for the mono-layer 
graphenei Two spatially inhomogeneous, spin singlet 
superconductors may arise in bilayer grpahene when elec- 
trons living on two layers attract each other. The order 
parameter is odd under the exchange of the layers, and 
breaks the translational symmetry of the lattice. Oth- 
erwise the s- and p-Kekule states are respectively, odd 
and even under the Dirac point exchange. At this mo- 
ment, the microscopic origin of the Nematic superconduc- 
tors are unknown. However, in a recent work, existence 
of some unconventional superconducting states has been 
proposed theoretically^ 



IX. SUMMARY AND DISCUSSION 

To summarize, we here present all the possible ordered 
phases, including fully gapped massive as well as ne- 
matic phases for fermions on bilayer graphene, and study 
their transformation under various symmetries (discrete 
and/or continuous). The parabolic bands touching each 
other at the Dirac points can be gapped out by forming 
either of the 28 massive ordered phases, which includes 8 



insulating and 4 superconducting phases. On the other 
hand, fermions can also condense into various nematic 
ground states. We here show that there as all together 56 
nematic order parameters in bilayer graphene. However, 
they corresponds to 16 semi-metallic and, 8 supercon- 
ducting nematic phases. Nematic superconductors can 
also break the translation symmetry. They can be clas- 
sified as nematic-Fulde-Farrell-Larkin-Ovchinikov super- 
conductors. Nevertheless, the realization of these cor- 
related ground states, the role of long ranged Coulomb 
interaction still remains as a open question. 

Here we have shown that in the low energy limit, the 
effective theory of gappless fermions can be described 
in terms of parabolic bands, touching each other at the 
Dirac points. However, upon taking into account the di- 
rect hopping amplitudes among the low-energy degrees 
of freedom (B sites), such parabolic band touching splits 
into four Dirac cones. Consequence, one still requires fi- 
nite strength of interactions to stabilize various ordered 
phases. Such critical strength of interactions is expected 
to be much smaller than that for monolayer graphene, 
since the direct hopping is weak. It is also therefore 
quite interesting to study the nature of the associated 
quantum phase transitions in bilayer graphene at zero 
temperature^ We leave this study for future investiga- 
tions. 
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XI. APPENDIX A: ALGEBRAIC DERIVATION 
OF MASSES AND NEMATIC ORDERS IN 
BILAYER GRAPHENE 

Here we present an alternate formulation to derive the 
number of massive and nematic order parameters in the 
bilayer graphene. In what follows next, our derivation 
closely mimics a similar analysis for monolayer graphene, 
presented in Appendix A of Ref. (53). However, the al- 
gebraic structures of these two problems enjoy significant 
differences, as we show below. Let us consider a 16 com- 
ponent Nambu-Dirac fermions "J = (^ p , ^ft) j so that 
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the effective Hamiltonian describing the low energy exci- 
tations reads as 



H k = H {k)®(-Hj(-k) 



where 



with 



Ho(k) = on dj 

i=l,2 



di = 



U2 _ u2 

2m 



and d-> = 



■'■ " y 



2m 



(41) 



(42) 



(43) 



oti and «2 are eight component Hermitian matrices.™ 
Here we wish to find all the Hermitian matrices (Mjs) 
that anti-commutes with Hamiltonian Hk, and develop 
gaps in the spectrum. In principle, there are numerous 
possibilities. However, we are are interested only in those 
matrices which give nonzero expectation values of the 
fermionic bi-linears 

m = ^ 0. (44) 

This condition can be satisfied only if 

M = - (<ti <g> J 8 ) M T ((7i <8 h) , (45) 

where Ig is eight dimensional unit matrix^ It was previ- 
ously shown by Altland and Zirnbauer— that there exist 
a unitary matrix U = U2 ® Ig, such that 

M = -M T , (46) 
where M = UMW. It particular one can show that 

U 2 = (±i) 1/2 e^V^V^-^ 3 , (47) 

whereas <j> = tt/4 was originally considered in Ref. (53). 
Therefore, after the unitary rotation by U, the mass- 
matrices are purely imaginary. Transformation of the 
free Hamiltonian Hq under the same unitary roatation 
can be captured by writing 

o,i = Re(cti) + i Im(a.i) . (48) 

Therefore 

H k = (CT 3 «'^e(a l ) + i a ®Im(a i ))d l = T l d l . (49) 

On the other hand, U^u^U^ = cr 2 . After the unitary 
transformation the free Hamiltonian is 

2 

H k = (02 ® Re(a z ) + i a <g> Im(oti)) dj = y^fjdj, 

i=l 

(50) 



with Ti now being purely imaginary matrices. A similar 
analysis showed that the free Dirac Hamiltonian (linear in 
momentum) in a monolayer graphene is defined in terms 
of real I\ matrices. 53 

Next we are after all the 16 x 16 imaginary matrices 
that anti-commute with 1^, i = 1,2. Since, iTj,iM are 
all real and square to — 1 , we first seek to know the max- 
imal number of q, so that for p > 0, the dimensionality 
of the real representation is sixteen and together close a 
C(p, q) algebra. The answer is 8. They form a Clifford 
algebra C(0, 8)i 55 ' 56 C(p, q) defines a set of p + q mutu- 
ally anti-commuting matrices; p of them square to +1, 
whereas q of them square to — 1 . Let us now define a set 
of 8 anti-commuting real matrices 



(51) 



all of which squares to — 1. Since, r 1; r 2 are imaginary, 
let us assume Ti = il\ and T2 — Hi- Next we want to 
find all the Hermitian imaginary matrices, anti-commute 
with I\ and Result is shown in the table below. 



Mass matrix(A'J) 



i h h{ = mi) 


1 


ilkhlmlnlp {k 7^ I ^ m ^ n ^ p = 3, • • • ,8) 


6 




6 


imi (I k IiI m I n ) (k ^ I ^ m ^ n = 3, • • • , 8) 


15 



The numbers in the right column indicates the number 
of mass matrices belongs to each class of mass matrices. 
Therefore we found there are all together 28 imaginary 
Hermitian matrices (mass orders) which anti-commute 
with the kinetic energy Hamiltonian. This number is in 
accordance with the one we compute here explicitly. One 
can perform a similar exercise to find the number of mass 
matrices for single layer graphene. That number is shown 
to be 3642 

One can also immediately find the number of nematic 
order in bilayer graphene. Any nematic order, with finite 
expectation value, assumes the form iYjM., where j = 1 
or 2. By construction it is Hermitian and imaginary. 
Therefore the total number of nematic order parameter 
is 2 (number of matrices in Hk) x28 (number of mass 
matrices) = 56, in agreement with our explicit computa- 
tion. 
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